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Abstract 

An 0(N) algorithm is proposed for calculating linear response func- 
tions of non-interacting electrons. This algorithm is simple and suit- 
able to parallel- and vector- computation. Since it avoids 0(N 3 ) com- 
putational effort of matrix diagonalization, it requires only 0(N) com- 
putational efforts where N is the dimension of the statevector. The 
use of this O(N) algorithm is very effective since otherwise we have to 
calculate large number of eigenstates, i.e., the occupied one-electron 
states up to the Fermi energy and the unoccupied states with higher 
energy. The advantage of this method compared to the Chebyshev 
polynomial method recently developed by Wang (L.W. Wang, Phys. 
Rev. B 49, 10154 (1994) ;L.W. Wang, Phys. Rev. Lett. 73, 1039 
(1994) ) is that our method can calculate linear response functions 
without any storage of huge statevectors on external storage. 



1 Introduction 

Computing the linear response functions (and the density of states) of large 
systems with thousands of atoms by using conventional methods requires 
to calculate the eigenvalues and eigenvectors of N x N Hamiltonian ma- 
trices (N ^> 10 6 ) from the lowest state to the Fermi energy and beyond 
it. The standard diagonalization routines are too much time-consuming in 
treating these problems because their computing time is proportional to N 3 
. Therefore efficient numerical algorithms, such as recursive Green's function 
methods |], ||, the Lanczos methods ||, J|, ^ ||, the Chebyshev polynomial 
expansion || [10], [□], [12|, [13], [14], |l|, [16|, [DJ, and conjugate gradient methods 
[0, §] have been developed and applied to various problems. 

In this paper, we present an efficient method for calculating the linear re- 
sponse functions of large quantum system. We give up the calculation of each 
exact eigenstates, instead we compute linear response functions by integrat- 
ing the time- dependent Schroedinger equations for a finite period determined 
by the required energy resolution. Since it avoids 0(N 3 ) computational ef- 
forts of matrix diagonalization, it requires only O(N) computational efforts 
for sparse Hamiltonian matrices. To realize this scheme, we exploit several 
numerical techniques such as the Chebyshev polynomial expansion of matrix 
functions i, 0, U, H (H, |]| ||, , random state vectors 0, H [jjj, 



Hamiltonian matrix discretized in real space [EG, 21 1, the time-dependent 
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Schroedinger equation discretized in real time [2!| ^3], [26|, |27], |28], p9 

ISO, BE B3, 031 PI BSf 



2 Time-dependent methods 

In this section, we describe how we reached the conclusion that we can cal- 
culate efficiently the linear response functions of large quantum systems by 
using the time-dependent homogeneous Schroedinger equations. 

2.1 Diagonalize or not diagonalize? 

Let us compare the computational efforts of the conventional diagonalization 
method and the time-dependent method by counting the number of floating 
point multiplications as a function of matrix dimension N, and show that 
time-dependent method is more efficient when large number of eigenstates 
are involved. 

First, we review the relation between the eigenstate representation and 
the time-dependent representation of linear response functions. The linear 
response function Xba^ + if}) of an observable B due to a monochromatic 
perturbation H ex = e -H w + M 7)* J 4 j s calculated by time-dependent perturbation 
theory f3T| . 

poo , ^ 

XBA {u + ir)) = {-i) dte +i ^ t {(E g \e +im Be- iHt A\E g )-c.c.} (1) 
« 2 f T dte +i{u)+iri)t \m[{E g \Be- iHt A\E g )e +iEst ) (2) 

J 

where \E g ) and E g are the groundstate of the many electron system and its 
energy, respectively; uj and r\ are the frequency and its resolution, respec- 
tively; T ^> 1/rj is the integration time. We use atomic units (a.u.), and 
indicate complex conjugate by "c.c". In numerical calculation of (Q), we 
have to discretize it in time, e.g., 

M 

Xba(u + i v ) = 2Y, Ate +i(u+iv)mAt lm[(E g \Be- iHmAt A\E g )e +iE9mAt } (3) 

m=0 

where M = T/At is the number of timesteps, T is the integration time in 
(0), and At is the width of timestep. On the other hand, we obtain the 
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eigenstate representation by inserting / = X)m=i \E m ){E m \ into 



(-5' g |-B|-£'-m)(-K m |A|-Eg) 
— ' ■ 1 i V ) - (E m - E g " 



loo 



m=l ^ 1 "'D ^g] 
_ A (-^g|^4|-^m)(-^m|-B|-5'g) 



(4) 



Next, we show the estimated computational efforts in Table [|. The diag- 
onalization method for N x N Hamiltonian matrix requires memory space 
of 0(N 2 ) and computational effort of 0(N 3 ). On the other hand, the time- 
dependent method requires memory space of 0(N 2 ) and computational effort 
of 0(MN 2 ) where M is the number of timesteps determined by the required 
energy resolution (See section |27|| ). By choosing an appropriate basis set, 
we can make the Hamiltonian a sparse matrix having only O(N) non-zero 
matrix elements 
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As the result, the computational effort and the 
memory space of the time-dependent method are reduced to O(MN) and 
O(N), respectively. Thus the time-dependent method can be more efficient 
than diagonalization method in large N limit. 



2.2 Newton or Schroedinger? 

Table |2| classifies various time-dependent methods in terms of kind of equation 
and homogeneity. Though the Newton equations of harmonic oscillators j37, 



39 



41] are mathematically equivalent to the Schroedinger equations 

1 



in the eigenstate representation, use of the Schroedinger equation |22|, 



30, 31. 32. 33. 34, 35] has advantage that we can exploit well 



developed concepts and formalism of quantum theory. It is true especially 
when we want to deal with quantum systems. Therefore, in this article, we 
study only the Schroedinger equations. 



2.3 Homogeneous or inhomogeneous? 

In this subsection we show that inhomogeneous time-dependent equations are 
more inefficient than homogeneous ones. This conclusion is valid not only 
for the Schroedinger equation (Particle Source Method p8| ) but also for the 
Newton equation ( Forced Oscillator Method J37], [38], [39], ||| £|l|) because 



both equations are equivalent in the eigenstate representation. 
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Let us define the computational effort of the time- dependent method by 
the number of timesteps M = T/At. Then the computational effort M 
is determined by the integration time T, because the maximum width of 
timestep is limited by the sampling theorem || independent of the detail of 
the method we use. The timestep should be much smaller than the inverse 
of band width, ir/E B , to reproduce the correct spectrum since, otherwise, 
according to @ we cannot distinguish the eigenvalues fl30 | 



2nk 

E k = E + — (* = 1,2,...). (5) 

In the following, we evaluate T for homogeneous and inhomogeneous 
Schroedinger equations to calculate the real-time Green's functions at many 
frequencies, Ui = lAu, I = 0, ±1,±2, within a required relative accu- 
racy 5. It turns out that T of inhomogeneous equations can be much longer 
than that of homogeneous equations. This conclusion applies to the linear 
response functions, too. 

First let us try to calculate the Green's function by solving the homoge- 
neous equation, 

l ±\frt) = H\frt) (6) 

with the initial condition |0; t — 0} = The auxiliary vectors are calculated 

as 

\kT) = H) F dt'\<j)-,t')e +i ^ +i W (7) 



o 



T 



H) / dt'e- iHt '\])e +l{wi+iri)t ' (8) 



o 

= (i-e^+^-^U) ( 9 ) 

« — — 7r\3) (10) 

UJl + 17] — H 

= Gfa + irilj) (11) 

where we have neglected the second term of (Q) by assuming T is large enough 
so that e~ r]T < 5. Therefore we estimate M for the homogeneous equation 

as 

T -log 5 



Mi 



At r/At 
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Next let us calculate the Green's function by solving the inhomogeneous 
Schroedinger equation, 

i±\frt) = H\cf>;t) + \j)(^]r e-^ + ^6{t) (12) 

with the initial condition \<f); t = 0) = 0. The solution at large T becomes 

|0; T) « Y,G(uj l + lV )\j)e-^' + ^ T (13) 
i 

where T satisfies e~ vT <C 5. Then the auxiliary vectors \<pi; T 2 ) are defined as 

1 & 



\4>V,T 2 ) = — f * dt'\frt)e-^' +i ^' 
±2 Jo 

= ^ Ptf^G{u* + iTj)\j)e-«'*- u ** f (14) 
T 2 Jo ■ 

= G(^ + ^)|j) 

i Ce -i ( w i -w i') T2 — l") 
+ £ G(^ + ^)|j) A— (15) 

« G(w,,+i77)|j> (16) 

where we have neglected the second term of fll5D by assuming that T 2 is large 
enough so that T 2 Auj 1/5. Therefore M becomes 

' (17) 



AcuAtS 

which can be much larger than Mi when Auj is small. 



3 Non-interacting Electrons 

In this section, we apply the time-dependent homogeneous Schroedinger 
equation to calculate efficiently the linear response functions and density of 
states of non-interacting electron systems, since it is well known that there 
exist wide and practically important areas in condensed matter physics where 
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non-interacting electron models are useful to predict various physical prop- 
erties. Hereafter we assume that the system is described by the one-electron 
Hamiltonian, 

H=-f + V(r). (18) 
3.1 Linear response function 

For non-interacting electrons, the linear response function (f|) can be rewrit- 



ten by using one-particle eigenstates as p6 



, v (i\B\j)(j\A\i} 
Xba{u) = 2^ 



E % <E^E 3 >E S ( w + ~ ( E 3 ~ E i) 

(t\A\j)(j\B\t) 

where Ef is the Fermi energy, and \i) and \j) are the occupied and empty 
one-particle states, respectively. This formula can be again rewritten in time- 
dependent representation as 

Xba{uj + if]) (20) 

= {-i) F dt J2 e +l(w+M?) ' {(i\e +iHt Be- iHt \j)(j\A\i) - c.c.) (21) 

Ei<Ef 
Ej>E f 

= (-i) ( T dt V e +i( - w+iv)t x 

{{%\9{E f - H)e + * m Be-* m 6(H - E f )\j){j\A\i) - c.c.} (22) 
' T dte +i{uj+iv)t K(t)\) (23) 



where the double brackets indicate the statistical average over random vectors 
|$), and K(t) is the time correlation function defined by 



K{t) = 2Im($|0(£ / - H)e +tHt Be- iHt 6(H - E f )A\<$>). (24) 



Equations fl23|) and (|24f ) are the main result of this paper. Note that calcu- 
lating the trace over the initial states \i) by using random vectors reduces 
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the computational effort by a factor of N. As the result, the computational 
effort still remains O(N) in spite of the double summation in flT9|) . 

In the above equations, we have introduced several numerical techniques. 
Firstly, the time-dependent statevectors, 

e~ im 9{H - E f )A\$) 
e- iHt 6{E f - H)\§) (25) 

are calculated by the leap frog method [^2], ^3], ^4], [25|, |2£ 

Ifat + At) = -2i At H\cj);t) + ](/>; t- At) (26) 

where the Hamiltonian matrix is discretized by finite difference p0| , |2l| 

pp. A N diff i 

H = E ^C^<f ) (x + nAx,y 7 z) + 0(Ax 2N ^f). (27) 

n=-N dlff 

Due to this discretization, the Hamiltonian matrix becomes sparse and the 
matrix vector multiplication in (|26|) can be done with O(N) computational 
complexity. We use N^ff = 4 formula in this paper. 

Secondly, the matrix step function for a normalized hermitian matrix X 
whose eigenvalues Xi are in the range [—1,1] is defined in its eigenstate basis 

6(X) = J2\X l ) 9{X i ) (Xi\. (28) 

By using this step function, we can avoid the difficulties in the partial sum 
in Operation of this function on an arbitrary vector |0) is numerically 



approximated by the Chebyshev polynomial expansion [0, 10, O, O, [L3L nl 



m y, 17 



f(X)\4>)*J2c k T k -i(X)\<t>) (29) 
fc=i 

where each term of the right hand side is calculated by vector recursion 
formulae 

T {X)\cf>) = |0) (30) 

T\(X)|0> = X\<f>) (31) 

T n+1 (X)\4>) = 2XT n {X)\ ( j ) )-T n ^{X)\ ( j>) n>l. (32) 



S 



To use this matrix function in (E3J), we should normalize the Hamiltonian 
matrix so that X = (H — Ef)/E norm has eigenvalues in the range [—1, 1]. 
Thirdly, we define random vectors with random phase by 

N 

|$) = \n)e +i4,n (33) 

n=l 

where \n) are basis vectors and —n < <p n < 7r, (n = 1, • • • , N) are uniform 
random variables that satisfy e -i<f> n i e iK yj _ Then we can derive 
various useful identities such as 

($|$) = £)<$|n)(n|$) =^e-^"e^« = JV (34) 

n n 

(( |$)($| )) = £>')(( e-^-'e*** )) (n| 

n'n 

= 5»<n| = I (35) 

n 



= £(n|A|n)=tr[A] = £(£j,4|£ m ) (36) 

n E m 

Equation ( |34] ) shows that each random vector is normalized to N, the number 
of one-particle eigenstates. Equation ([35|) shows that random vectors have 
normalized completeness. Equation (|36| ) shows that the expectation value 
of an operator by random vectors gives the trace of the operator. We used 
this identity to calculate the trace over \i) in (^3|) and ([Zip. These random 
vectors with random phase are more useful in calculating expectation values 
than random vectors with random amplitude since each random vectors are 
automatically normalized. 

Finally the formula for numerical calculation of polarizability function 
Xpa{u) with a, (3 = x,y, z becomes 

X/*aM « [i^ll dte-^{e +i ^-8 Pa )K{t)\\ (37) 

K(t) = v{ jl ir]) M mEf - H)e+ iHt x 

p^ iHt d{E mt - H)9{H - E f )p a \$) (38) 
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where V is the volume of the supercell, and the dipole moment operators 



(j\A\i) = (j\x a \i) (39) 
(i\B\j) = ^(iM), (40) 

are modified to momentum operators by partial integration. We also inserted 
a low energy projection operator 6{E cut —H) into ( |38] ) to eliminate unphysical 
high energy components of the random vectors. This filter is much more 
effective than the quadratic filter used in |T3. In calculating very large 



systems, we need only few random vectors for statistical averaging, since the 
fluctuation becomes smaller as the system size N becomes larger ||28|| . 

Figure p] shows the dielectric function e xx {oj) = 1 + Aixxxxi^) of four 
electrons in three dimensional harmonic potential 

V(r) = ^ (41) 

calculated with 32 3 cubic meshes, ujq = 0.1, r\ = 10~ 4 . Three random vectors 
are used. The analytical result JE| 



e xx (u) = l + —A— 2 ^— 42 

V UJq — CO 2 — IUT] 

is also shown for comparison, where N e is the number of electrons in the 
supercell of volume V. The deviation from the exact result near uo = is due 
to finite rj. The result shows that our method works very well for uj 3> r]. 

Figure [| shows the dielectric function with energy resolution rj = 0.05(eV^) 
of silicon crystal consisting of 2 15 Si atoms in a cubic supercell of 16 3 unit cells. 
Each unit cell is divided into 8 3 cubic meshes. One random vector is used. 



We used the empirical local pseudopotential in reference i43| . The result 
agrees with experimental results and other theoretical calculations [13, |3Jj . 

In some cases we may want to ask which part of the real space the elec- 
trons contributing to the linear response function come from. We can answer 
to this question by calculating the linear response function by restricting the 
range of the trace in (|3"T|) within a real space domain D. This can be done 



by replacing |$) by |<3>') = Pd|<&) where Pjj = J2 n <=D \ n )i n \ 1S the real space 
projection operator onto D. 
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3.2 Density of states 

The density of states of the system can be calculated as |32 

p{uj) = — V Im G nn {u + irj) = — Im (tr [G(u + irj)}) (43) 

7T V 7T 



by combining (|TT| ) and (p6|). 

Figure [5] shows the numerical and analytical results of the density of states 
in 3D harmonic potential with 32 3 cubic meshes, ujq = 0.1, and rj = 10~ 3 . 
Three random vectors are used. Figure |] shows the density of states of 
silicon crystal consisting of 2 15 Si atoms in a cubic supercell of 16 3 unit cells. 
Each unit cell is divided into 8 3 cubic meshes. The energy resolution is 
7] = 0.05(eV). Three random vectors are used. 

We can also calculate the local density of states integrated in a given 
domain D by using the real space projection operator Pp to restrict the 
summation in (|33[) within D, 



Pd{u) = — E Im G nn (u + ir}) = — Im (\x[P D G(u +irj)\) . (44) 



7T r- n 



Photonic band structures in two-dimensional periodic structure of dielec- 



tric material j|6], eH can a ^ so ^ e calculated by using (fE|) or (|44]) since the 
Maxwell equations of this system are reduced to the Schroedinger equation 
with position dependent mass, i.e. , 

HH z (x,y) = ^H z (x,y)=EH z (x,y) (45) 

d -1 d d -1 d 

dx e(x, y) dx dy e(x, y) dy 

for if-mode where H z is the z-component of the magnetic field, and 

HE z {x,y) = —E z (x,y) = EE z (x,y) (47) 




^ e(x,y) \ dx 2 dy 2 ] 

for S-mode where E z is the z-component of the electric field. 

Figure [5] shows a typical structure of two-dimensional photonic crystal 
cavities used in our calculation, and Fig. |6] shows the calculated density of 
states as a function of frequency and wave number. 
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4 Summary 



In this article we proposed a new numerical method suitable for calculating 
the linear response functions (and the density of states) of non-interacting 
electrons, in which the sum over the initial one-particle states are efficiently 
calculated by using random vectors. The advantage of this method compared 
to the Chebyshev polynomial method by Wang to calculate optical absorption 
of non-interacting electrons [14j] is that our method can calculate not only the 
imaginary part but also the real part of the linear response functions at the 
same time, and that it can calculate them without any input-output of stat- 
evectors on external storage. As the result, our method can calculate much 
larger systems than Wang's method. The Chebyshev polynomial method of 
degree M should store O(M) statevectors of size O(N) on external storage 
to make the table of 0(M 2 ) generalized Chebyshev moments A m>m i and may 
take very long I/O time. 

The application of this method to photonic band structures, silicon nanocrys- 
talites, and periodic structures of chaotic systems will be presented elsewhere 

EM, El- 
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Diagonalization Method 



^ (Eg\B\E m )(E m \A\E g ) " (E g \A\E m ) (E m \B\E g ) 
ti fa + iv) ~ (Em - E g ) ^ (u + irj) + (E m - E g ) 



Calculation 


Dense Matrix 


Sparse Matrix 


Computation 


Memory 


Computation 


Memory 


E \E ) 

J -'mi \ J - J ml 


TV 3 


N 2 


N :i 


N 2 


(E 9 \A\E n ) 


N 2 


N 


N 


N 


E 

m 


N 


N 


N 


N 



Time-dependent Method 

M 

XBA(w + ir]) = 2 &te +t{uJ+E ° +tri)mAt Im(E g \Be- iHmAt A\Eg) 

m=0 



Calculation 


Dense Matrix 


Sparse Matrix 


Computation 


Memory 


Computation 


Memory 


e- lHmAt B\E g ) 
(Eg\Ae~ iHmAt B\E g ) 

M 

E 

m=0 


MN 2 
N 2 

M 


N 2 
N 

1 


MN 
N 

M 


N 
N 

1 



Table 1: Comparison of diagonalization method and time-dependent method 



17 





equation 


homogeneous 


inhomogeneous 


classical mechanics 


Newton 


harmonic oscillator 


forced oscillator 
[0 0, 0, © 


quantum mechanics 


Schroedinger 


r 


IDSE 

22, 23, 24, 25 

n, m, n, n 


MM 


particle source 
|28| 



Table 2: Comparison of time-dependent equations 
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Figure 1: e xx (uj) of four electrons in a three dimensional harmonic oscillator 
calculated with 32 3 cubic meshes, ujq = 0.1, 77 = 1CT 4 ; (a) real part, (b) 
imaginary part. 



Figure 2: e xx (u) of silicon crystal consisting of 2 15 Si atoms in a cubic su- 
percell of 16 3 unit cells. Each unit cell is divided into 8 3 cubic meshes. The 
energy resolution is r\ = 0.05(eV). We used the empirical local pseudopoten- 
tial in reference (a) real part, (b) imaginary part. 



Figure 3: Density of states of 3D harmonic oscillator calculated with 32 3 
cubic meshes, u = 0.1, and r\ = 10~ 4 , and analytical result. 



Figure 4: Density of states of silicon crystal consisting of 2 15 Si atoms in 
16 3 unit cells. Each unit cell is divided into 8 3 cubic meshes. The energy 
resolution is 77 = 0.05(6^). 



Figure 5: A typical structure of two-dimensional photonic crystal cavities 
used in our calculation 



Figure 6: The calculated density of states as a function of frequency and 
wave number 
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